# -*- coding: utf-8 -*-
"""
    @project: pythonProject
    @Author：HanYonghua
    @file： Decimation.py
    @date：2025/4/22 17:26
    @blogs: https://www.ncatest.com.cn
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# 模拟250 MSPS采样率的信号
fs_original = 250e6  # 原始采样率 (250 MHz)
t = np.arange(0, 1e-5, 1/fs_original)  # 10μs数据
f_signal = 10e6  # 信号频率 (10 MHz)
adc_data = np.sin(2 * np.pi * f_signal * t) + 0.1 * np.random.randn(len(t))  # 含噪声

D = 5  # 抽取因子，目标采样率 = 250 MHz / 5 = 50 MHz
cutoff = fs_original / (2 * D)  # 截止频率 = 25 MHz (Nyquist频率)
taps = signal.remez(101, [0, 0.9 * cutoff, 1.1 * cutoff, fs_original/2], [1, 0], fs=fs_original)
filtered_data = signal.lfilter(taps, 1.0, adc_data)
 
decimated_data = filtered_data[::D]  # 每D个样本取一个
fs_new = fs_original / D  # 新采样率 (50 MHz)

# 绘制频谱对比
def plot_spectrum(data, fs, title):
    fft_data = np.fft.fft(data)
    freqs = np.fft.fftfreq(len(data), 1/fs)
    plt.plot(np.fft.fftshift(freqs)/1e6, 10*np.log10(np.abs(np.fft.fftshift(fft_data))**2))
    plt.title(title)
    plt.xlabel("Frequency (MHz)")
    plt.ylabel("Power (dB)")

plt.figure(figsize=(12, 8))
plt.subplot(2, 1, 1)
plot_spectrum(adc_data, fs_original, "Original Signal (250 MSPS)")
plt.xlim([-fs_original/2e6, fs_original/2e6])
  
plt.subplot(2, 1, 2)
plot_spectrum(decimated_data, fs_new, f"Decimated Signal ({fs_new/1e6} MSPS)")
plt.xlim([-fs_new/2e6, fs_new/2e6])
plt.tight_layout()
plt.show()
